Universality of the melting curves for a wide range of interaction potentials 
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We demonstrate that the melting curves of various model systems of interacting particles collapse 
to (or are located very close to) a universal master curve on a plane of appropriately chosen scaled 
variables. The physics behind this universality is discussed. An equation for the emerging "universal 
melting curve" is proposed. The obtained results can be used to approximately predict melting of 
various substances in a wide range of conditions. 

PACS numbers: 05.70.Fh, 64.60.-i, 64.70. D-, 64.70.dj 



Exploring phase behavior of different substances is an 
outstanding physical problem with significant impact on 
basic (e.g. theories of phase transitions) and applied 
(e.g. materials science) research. Significant progress in 
this area has been achieved over the last decades using a 
variety of theoretical, experimental, and computational 
methods. This resulted in detailed phase diagrams of 
some conventional atomic and molecular substances, as 
well as different kinds of new materials such as fullerenes 
and soft matter (colloidal suspensions, polymers, surfac- 
tants, complex plasmas, etc.) [H-Q- 

The traditional procedure of predicting the phase dia- 
gram of a substance is based on extensive numerical sim- 
ulations using a model pair potential which approximates 
actual interactions in this substance. Simple pair inter- 
actions studied in this context include the hard sphere, 
inverse power law, large family of Lennard- Jones-t ype , 
Yukawa, Gaussian and other model potentials |5l4ll|. 



There exists also several approximate methods to lo- 
cate phase boundaries (especially for the fluid-solid phase 
changes), which are independent of an exact shape of 
the interaction. These are well known phenomenolog- 
ical criteria for freezing and melting, like e.g. Linde- 
mann melting law, Hansen- Verlet freezing rule, Raveche- 
Mountain-Street criterion for freezing, and a dynamical 
criterion for freezing in colloidal suspensions (for a re- 
view see Ref. |12|). These criteria are typically based on 
the properties of only one of the two coexisting phases 
and predict quasi-universal values of certain structural 
or dynamical quantities. Quasi- universality means that 
a quantity is not exactly constant, but varies in a suffi- 
ciently narrow range for a broad variety of physical sys- 
tems (i.e. interaction shapes). These useful empirical 
rules can sometimes be directly applied to the substance 
under investigation. More often, they are used in com- 
bination with numerical simulations, which allows to ap- 
proximately locate phase boundaries with a very modest 
computational cost [l^. Nevertheless, it would be often 
desirable to have simple analytical expressions describ- 
ing phase coexistence, based only on the properties of 
the interaction potential. Even though such expressions 
can be not completely universal (e.g. their applicability 
limited to a certain class of interactions), they can be 



quite helpful, especially in cases when interparticle inter- 
actions depend on a number of system parameters, which 
can vary from one situation to another [3]. 

The main purpose of this Letter is to describe a univer- 
sality of melting curves for a wide range of interparticle 
interactions. In particular, we identify a pair of scaled 
variables, which characterize the system of strongly in- 
teracting particles. We then show that in a plane of this 
variables the available numerical data for melting of quite 
different model systems are located on (or very close to) 
a single master curve. A simple expression for this "uni- 
versal melting curve" is proposed. We briefly discuss this 
finding in the context of several relevant recent studies. 

The physical idea is rather simple and transparent. Let 
us consider a system of particles interacting via a pairwise 
interaction potential exhibiting conventional properties: 
The interaction force is strongly repulsive at short sep- 
arations and vanishes at infinite separations (the long- 
range asymptote can be either repulsive or attractive). 
We then assume that the particle system is strongly cou- 
pled, i.e. the system is either in a dense fluid or solid 
state. The particles form a regular structure where large 
deviations of the interparticle separation from its average 
value are very seldom. Consequently, the state of the sys- 
tem should be virtually insensitive to the exact shape of 
the interaction potential at short distances. If, in addi- 
tion, the potential decays sufficiently fast for distances 
beyond the mean interparticle separation, the proper- 
ties of the system are hardly influenced by the shape of 
the interaction potential at long distances. The behav- 
ior of the potential at the average interparticle separa- 
tion plays a dominant role. Physically, the force and its 
derivative evaluated at the mean interparticle distance 
are the two quantities which should really matter. We 
therefore put forward the following hypothesis: Two sys- 
tems of strongly interacting particles, characterized by 
different interaction potentials Ui{r) and U2{r) should 
exhibit similar properties provided U[{A) — U2{^) and 
C/"(A) ~ J/j'l^); where A is the mean interparticle dis- 
tance. The purpose of this Letter is not to give a general 
proof of this hypothesis, but demonstrate that it leads 
to remarkably good results when applied to fluid-solid 
coexistence of various strongly coupled systems. 
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To proceed further it is useful to chose a reference po- 
tential of interaction. Since fluid-solid phase changes are 
dominated by the repulsive portions of the interparti- 
cle forces, the natural and simplest choice is the inverse- 
power-law (IPL) family of potentials U{r) = £(cr/r)", 
where e and cr are energy and length scales. This po- 
tential defines a system of "soft" spheres and is used to 
model e.g. one-component plasma (OCP) and simple 
metals under extreme thermodynamic conditions. The 
softness of the potential is governed by the index n and 
can be smoothly varied in essentially entire range: the 
very soft n — I repulsion corresponds to the OCP limit, 
the potential becomes harder as n increases, and for 
n oo it approaches the hard sphere (HS) limit. 

In the absence of external fields and dissipation the 
system of N particles interacting via the IPL potential 
is described by Newton's equations of motion, which in 
reduced units have the form 
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(1) 



where T — rt(e/T)(cr/A)" = nT^^p'l^^ characterizes the 
force of interaction at the mean interparticle distance, 
A = iy I'N'f'^'^ is the structure-independent interparticle 
spacing, V is the system volume, — T/e is the reduced 
temperature and = Na'^ /V is the reduced density. 
The single scaled variable J- [or any other combination 
of /5*^^/T* with/without n] describes the whole range of 
thermodynamic properties of the given IPL potential. 

The phase diagram of the IPL system has been exten- 
sively studied 0,11, 15, Hi- Three phases can exist, one 



fluid and two solid forming either a body-centered cubic 
(bcc) lattice, favorable for soft interactions < 6), or a 
face-centered cubic (fee) lattice, favorable for hard inter- 
actions (n > 6) 0, Q • Since we are mostly interested in 
the fluid-solid transition here, the phase diagram of the 
IPL system in the (1/n, J") plane plotted in Fig.[T]shows 
only the numerical data for freezing and melting. 

The data in Fig. [1] are from numerical simulations of 
Refs. 0, [lE [3i and cover the wide range of softness 
between n ~ 3 and n = 100. Note that when n ap- 
proach 3 from above (in three-dimensional case under 
consideration), the difference between freezing and melt- 
ing densities vanishes and melting and freezing curves 
fall on top of each other. When n < 3 the thermody- 
namic functions of the IPL system are undefined, unless 
a uniform "neutralizing background" is provided. In this 
case the fluid-solid phase change occurs at constant den- 
sity (not at constant pressure), and freezing and melting 
points would coincide in the plane (1/n, F). Two rele- 
vant simulation points, corresponding to the OCP regime 
at n = 1 are also shown in Fig. [TJ 

Having defined and discussed the properties of the ref- 
erence potential, let us identify the convincing pair of 
scaled variables, which will be used in further analy- 
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FIG. 1: (Color) Phase diagram of the IPL system in the T 
vs 1/n plane, where T = n{e/T)(a / A)"^ is the interaction pa- 
rameter. Here red (blue) squares, circles, and triangles corre- 
spond to the representative numerical data for freezing (melt- 
ing) taken from Refs. Q, and [l^, respectively. Diamond 
and star mark the OCP in = 1) fluid-solid phase change as es- 
timated in Ref. Q] and [l3|, respectively. Red (blue) dashed 
curves correspond to the hard-sphere asymptotes (n oo) 
for freezing (melting), see text. 



sis. We require that the first and the second deriva- 
tive of some arbitrary potential U{r) and the IPL po- 
tential are equal at r = A. This immediately yields 
n = -1 - U"{A)A/U'{A), which is one of the possi- 
ble definitions of a local effective IPL exponent [l^ . For 
a given n the phase state of the IPL system is fully deter- 
mined by the value of the parameter T = —U'{A)A/T. 
This provides us with the required pair of scaled vari- 
ables: the generalized softness parameter 
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and the generalized interaction (force) parameter 

T ^ -U'{A)A/T. (3) 

Our hypothesis can now be slightly reformulated in the 
form of the approximate corresponding states principle: 
Two different systems of strongly interacting particles 
having the same values of s and J^behave alike. An ideal 
way to verify the vitality of this principle is to study the 
fluid-solid coexistence of various physical systems. 

We consider a wide range of interactions including 
purely repulsive and those having long-range attractive 
branches. It is convenient to present the potentials 
in the general form U{r) = eu(x), where x = r/a 
is the normalized distance. In this notation we have 
u{x) — x~" for the IPL potentials. Another purely re- 
pulsive potential is the Yukawa (Debye-Hiickel) potential 
u{x) — (l/x) exp{~x) extensively used in the context of 
colloidal suspensions, conventional electron-ion plasmas, 
and complex (dusty) plasmas [11-11]. One of the most 
widely studied model potential capturing much of the es- 
sential physics of simple atomic substances is the 12 — 6 



Lennard-Jones (LJ) potential, u{x) = 4(a; 
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FIG. 2: (Color) Representative set of pair interaction poten- 
tials considered in this Letter. See text for description. 



There are also several modifications to the conventional 
12 — 6 LJ potential, including n — 6 family of potentials, 

u{x) = - 6)](n/6)S/("-^)(a;-" - x"^), 

where 6 < n < 12. The exp —6 potential is defined as 



u(x) 



-oo. 
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where x-ca corresponds to the maximum of the function 
in the second line, the repulsion is of "hard sphere" type 
for X < Xm and is controlled by the coefficient a for 
X > Xm- The Weeks-Chandler- Anderson (WCA) refer- 
ence repulsive potential correctly describes the repulsive 
forces in the 12 — 6 LJ fluid for x < 2^/^ and is set equal 
to zero at longer distances One more potential we 

consider has a Gaussian form, u{x) = exp(— x^) and is 
referred to as the Gaussian core model (GCM) [20j . It is 
widely used to describe effective interactions in soft mat- 
ter physics (e.g polymers in solutions) Note that the 
bounded GCM potential does not diverge at a; as do 
the other potentials considered here. We have plotted a 
representative set of pair interaction potentials in Fig. [2] 
Figure [3] summarizes numerical data related to melt- 
ing of the model potentials discussed above in the (s, 
J-) plane. The data for the IPL potential are the same 
as in Fig. pQ data for the Yukawa potential are from 
Refs. [i, HJll^l , for the 12 — 6 LJ potential we use nu- 
merical data from Refs. (gI. [23l - l26j | . for the n — 6 potential 
we use data from [2^, exp— 6 data for a = 13 are from 
23, WCA data are fromJilLand the data for the GCM 
potential are from Refs. [16|,l29|. Amazingly, most of the 
data points fall on (or very close to) a single curve - the 
"universal melting curve" . No systematic deviations are 
evident, except for the case of GCM potential (although 
even in this case data points are still located in the prox- 
imity of the universal curve) . This is not surprising since 
GCM has a soft repulsive core and, therefore, for high 
densities one may expect some dependence of the system 
properties on the exact shape of the potential at short 
distances (absent for diverging repulsion). 
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FIG. 3: (Golor) The "universal melting curve" in the (s, T) 
plane. Symbols (defined in the figure) denote the available 
numerical data. The solid curve is the fit of Eq. Q. 



We can propose a simple fit for the emerging "univer- 
sal melting curve" . For sufficiently soft interactions [soft 
spheres (SS) with 0.1 < s < 1] the data points can be 
reasonably described by a simple power-law dependence 
J"ss — 106s^/'^ (dashed line in Fig. where the front 
factor 106 is taken to agree with the OCP (s = 1) value 
from Ref. Q. In the opposite limit of HS-like interac- 
tion (s — s- 0) we use the following arguments (applied 
previously to GCM, e.g. (l^l)- The effective hard-sphere 
diameter a;* can pe approximated as a distance at which 
the pair-interaction Boltzmann factor exp[— eM(a:)/T] is 
equal to 1/2, i.e. u{x^) = TH,ln2. The HS volume frac- 
tion at melting is {n/6)xlp^ ~ 0.545 (~ 0.494 at freez- 
ing). For the IPL system T = np"^^/T^,, which imme- 
diately yields T ~ n(1.041)"/'^ ln2 at melting [similarly 
J" ~ n(0.943)"/^ln2 at freezing]. This asymptotic esti- 
mates for melting and freezing are shown in Fig. [T] by 
the dotted curves. They approach numerical data at low 
softness (s < 0.02). Generalizing the IPL case to other 
interactions we get the following HS asymptote for the 
melting curve Tns ^ (ln2/s)(1.041)l/3^ The simplest 
form which is correct in the corresponding limiting cases 
and yields reasonable interpolation between them is 



^=(J-S^S + -^HS) 



(4) 



Equation ^ with i' — 8/5 do a remarkably good job (see 
Fig. El . 

Let us briefly discuss the observed universality in the 
context of some previous studies related to this topic. 
It was reported in Ref. |1J| that the freezing curves of 
Yukawa and IPL potentials exhibit striking similarity 
when plotted in the plane of appropriately chosen param- 
eters (the same would apply for the melting curves, espe- 
cially when the interactions are sufficiently soft, so that 
freezing and melting curves are almost indistinguishable) . 
The pair of generic parameters identified in Ref. [141] is 
the actual coupling strength [T = U ( A)/T] and the steep- 
ness of the interaction potential at the mean interparticle 
distance [7 = |c? In C/(r)/(i In r|r=A]- This choice comes 
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from requiring the potential itself and its first derivative 
to be equal to those of the IPL potential at r = A, al- 
though this was not stated explicitly in ■ Our present 
approach based on the equality of the force and its first 
derivative is more physically sound (for instance, it al- 
lows us to deal with attractive interactions) and is thus 
superior to that of [31 ■ However, it is easy to see that for 
Yukawa and IPL potentials both approaches yield similar 
results. Let us require that T and 7 are equal for Yukawa 
and IPL potentials. Then T are equal too, but the gen- 
eralized softness parameters s are different by the factor 
n^/ (n^ — n + 1), which is 1 for n = 1, tends to 1 for large 
n, and reaches the maximum value of ~ 1.3 at n = 2. 
This maximum discrepancy between the two approaches 
corresponds to the regime where no numerical data for 
the IPL potential are presently available (3 > n > 1). 

Accurate approximations for the dependence of tem- 
perature on density along the freezing and melting curves 
of the 12-6 LJ fluid have been recently proposed (30| . 
Mehing equation has the form ~ 1.97 pt - l.OSpl. It 
is essentially "exact" in the high temperature limit and 
at the triple point and shows remarkably good agreement 
with the numerical data in the intermediate region. Sim- 
ilar equation can be easily derived on the basis of the 
present consideration. Note that the softness of the 12-6 
LJ potential varies in a relatively narrow range between 
~ 0.083 in the high-temperature limit and ~ 0.053 at 
the triple point. In this range the universal melting curve 
exhibits minimum and the dependence of J-" on s is rather 
weak. In the first approximation we can simply assume 

~ const, where const ~ 24 (see Fig. [3]). For the 12-6 
LJ potential this yields ~ 2.0p^ — l.Op^, in reasonable 
agreement with a more involved consideration [30| . 

Finally, we would like to point out that the reported 
universality of melting curves can be considered in the 
context of the so-called "isomorph" concept (see Ref. [3l| 
and references therein). In particular, we have seen that 
melting of a wide class of strongly coupled systems can be 
well reproduced by a properly chosen IPL reference sys- 
tem. More generally, our results are not inconsistent with 
the prediction that there exists a class of fluids which ex- 
hibit equivalence principle (in the sense that the system 
of isomorphic curves of one fluid can be mapped onto 
that of the other fluid) and that the fluid-solid coexis- 
tence curve is itself an isomorph [31]. 

To summarize, we have shown that the melting curves 
for a wide range of different potentials exhibit a universal 
shape on a plane of appropriately chosen reduced vari- 
ables. These are the normalized force of interparticle 
interaction and the generalized softness parameter, both 
evaluated at the mean interparticle distance. Available 
numerical data for the melting of IPL, Yukawa, WCA, 
12-6 LJ, n — 6 LJ, exp —6 systems are essentially collaps- 
ing on a single curve, data for the bounded GCM poten- 
tial are not exactly on the curve, but in its proximity. 
We propose an analytical expression for this "universal 



melting curve". Although this phenomenological melt- 
ing equation cannot replace the proper thermodynamic 
prescription of the phase portraits of different systems, 
we believe that it can be very useful in approximately 
predicting, with very little effort, the melting transitions 
for various substances over a wide range of conditions. 
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